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SPECTRAL  CHARACTERISTICS  OF  A FLICKER  NOISE  MODEL 


1 . INTRODUCTION 


Various  papers  have  been  written  on  the  simulation  of  "flicker"  or  "y" 
noise.  However,  the  majority  of  them  are  concerned  with  the  effect  of 
flicker  noise  on  frequency  measurement  and  hence  Investigated  flicker  noise 
as  a phase  variation.  We  are  more  Interested  In  simulating  flicker  noise 
Inherent  In  the  output  signals  from  Infrared  detectors  and  their  associated 
electronics. 

There  are  several  possible  methods  of  slmllatlng  flicker  noise.  One 
method  Is  to  consider  a filter  operating  on  white  noise  which  gives  the 
required  spectral  variation.  This  Is  basically  the  method  Barnes  and 
Allan  (1)  used  In  obtaining  an  expression  for  variation  In  phase  due  to 
flicker  noise.  Another  system  which  would  seem  to  have  a more  reasonable 
basis  as  a physical  model  for  flicker  noise  Is  the  mechanical  model  of 
Halford  (2)  which  considers  classes  of  time-dependent  perturbations  occurring 
at  random,  generating  random  noise  with  the  required  spectral  density. 


2 . THEORY 


We  elected  to  use  the  simpler  approach  of  a filter  operating  on  white 
noise  as  shown  In  Figure  1. 

We  wish  to  generate  noise  with  the  spectral  power  density,  <|)(f),  of  the 

form 


■Kf) 


. (1) 


where  X = for  "flicker"  noise. 


If  we  pass  white  noise,  with  a constant  spectral 
a network  with  a power  transfer  function  proportional 


power  density,  through 
to  — ^ , then  the 

f^A 


1 


output  spectrum  will  have  the  required  characteristics.  The  network  must 
have  a signal  transfer  function  of  the  form 


where  K is  a constant.  However,  a digital  computer  simulation  will 
generally  be  working  in  the,  time,  rather  than  the  frequency  domain,  so  it 
is  necessary  to  define  the  processing  network  G(f)  in  the  time  domain. 

The  impulse  response  H(t)  of  C(f)  is  its  Fourier  Transform,  that  is 


00 

-f 


G(f)e'^  ^’^df 


It  is  shown  in  Appendix  I that  a possible  form  is 


2K  sin(rX)  r(l-X)(2irt)  , t>o  and  o<X<l  ...  (4) 


/-X  X-1 

ex  dx  for  X>0.  The  output  signal  y(t)  in  response 


to  x(t)  is  normally  given  by  the  convolution  Integral 


y(t)  = / x(t)  H(t-x)  dt 


...  (5) 


However,  this  Integral  does  not  exist  when  x(t)  is  white  noise  and  the 
Impulse  H(t)  is  of  the  form  given  in  equation  (4).  To  overcome  this 
problem  we  will  modify  the  Impulse  response  so  that  It  Is  non-zero  for  only 
a finite  duration  T and  is  zero  outside  this  range.  Such  a change  amounts 
to  modifying  the  power  transfer  function  so  that  it  does  not  follow  an 
_2X 

f law  down  to  zero  frequency,  but  introduces  a levelling  off  at  low 
frequencies.  The  output  y(t)  will  now  be  expected  to  have  the  required 
spectral  law  only  over  a restricted  frequency  range  and  as  such  can  be  best 
described  as  a quasi-flicker  noise  process.  The  modified  (truncated) 
Impulse  will  be  taken  as 


H^(t)  = 2K  sin(TtX)  F (1-X)  (2Trt) o<t$T 
H^(t)  = 0 elsewhere 


t 


Equations  (5)  and  (6)  are  applicable  to  continuous  processes.  However,  for 
digital  computer  simulation  a discrete  process  with  sampling  at  regular  time 
Intervals  At  must  be  used.  Therefore  we  consider,  Instead  of  the  con- 
tinuous Hj^(t),  a discrete  value  Impulse  response  which  at  time  jAt  Is 


j 


Hj^(jAt)  - h^6(t  - lAt)  , 


1»1 


J^or  J«N  where  NAt  “ T,  the  duration  of  the  Impulse  response. 
(jAt)  Is  zero  as  required  by  equation  (6). 

Using  (7),  and  a discrete  form  of  (5),  It  follows  that 


j-1 

I'j  “ y(jAt)  - At  ^i^j-1 

1-j-N 


If  we  now  consider  the  particular  case  of  A = *s,  then 


...  (7) 


For  j>N, 


...  (8) 


H^(t)  - VI  Kt“** 


...  (9) 


The  coefficients  h^  Introduced  In  (7),  can  be  defined  in  two  ways.  Firstly 
It  can  be  defined  that 


hj^At 


(l+l)At 


...  (10) 


lAt 


for  1 Integer  in  the  range  0 i < N.  This  approach  leads  to 


j = 2vT  K(At)“^  ^1  + 1 - yfi^  ...  (11) 


Rather  than  the  previous  approach,  we  choose  to  define 


= /T  K(lAt)‘ 


..  (12) 
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Equation  (8)  is  evaluated  using  the  N coefficients  h,  to  h„  defined  by  (12) 
to  give 


I 


...  (13) 


This  equation  then  uses  only  the  last  N values  of  x to  generate  y.(N). 

^1  ^ 

This  then  introduces  a low  frequency  cutoff  of  the  j behaviour  at  a 

1 

frequency  of  about  . The  upper  frequency  cutoff  is  at  the  Nyquist 
frequency  of  . 


The  problem  now  arises  of  investigating  how  closely  the  expression  in 

equation  (13)  produces  a spectral  power  density  with  a ^ characteristic  over 

the  expected  band.  There  are  at  least  two  methods  of  estimating  spectral 
power  density. 


Firstly,  an  estimate  of  the  power  spectrum  can  be  obtained  via  the 

Fast  Fourier  Transform  (3)  (FFT) . However,  as  stated  by  Furber  (4)  a 

single  FFT  calculation  from  one  set  of  samples  of  a random  signal  is  a poor 

estimate  of  the  power  spectrum  and  is  a statistically  unstable  estimate. 

The  stability  can  be  increased  by  performing  FFT  calculations  on  a number 

(M)  of  different  sets  of  samples  of  the  stationary  signal  and^^then  average 

these  together.  The  variance  of  the  estimate  decreases  as  — . 

M 

The  second  method  is  to  calculate  the  autocorrelation  function  and  then 
calculate  the  Fourier  transform  of  this  function.  Since  we  could  calculate 
an  analytical  solution  for  the  autocorrelation  function  we  used  this  second 
method  for  calculating  the  spectral  power  density. 


3.  CALCULATION  OF  AUTOCORRELATION  FUNCTION 


It  is  shown  by  Papoulis  (5)  that  the  autocorrelation  function  of  the 
random  signal  y(t)  generated  by  applying  white  noise  of  uniform  spectral 
power  density  S to  a system  with  impulse  response,  H(t),  is  given  by 


k\ 

V 


00 

y*  H(u)  H( 


R(t)  = S J H(u)  H(t+u)  du  , 
o 


...  (14) 


where,  as  usual. 


R(t)  = < y(t)  y(t-T)  > , 


...  (15) 


I 


where  t is  the  lag  variable  and  the  brackets  signify  expectation  value. 

It  is  readily  shown  that  the  discrete  form  of  (14)  can,  for  a finite 
duration  impulse  response,  be  written  as 

N-n 

R(n,N)  = SAt  h,  h . . , ...  (16) 

j n + 3 
j=l 

Where  the  delay  time  t is  given  by  nAt.  It  should  be  noted  that  the  above 
expression  was  obtained  for  0 < n < N - 1.  However,  it  is  readily 
demonstrated  that 


I 


R(n,N)  = R(|n|,N) 


for  - (N-1)  n 0. 

These  latter  results  can  be  combined  to  give 


(17) 


R(n,N) 


SAt 


N-|n| 

L 

j=l 


+ j 


(18) 


for  - (N-1)  < n ^ (N-1). 

Using  the  h coefficients  already  defined  for  X = yields 


N-|n|  ^ 

R(n,N)  = 2 SK^  ■- 

^ /j(j  + |n|) 


or  expressed  differently 


R(n,N) 


2 SKV(n,N)  , 


N-|n|  ^ 

where  R^(n,N)  = N — 

+ |n|) 


I 

i 


...  (19) 


J 
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4.  CALCULATION  OF  SPECTRAL  POWER  DENSITY 


For  a continuous  process,  the  spectral  power  density  corresponding 

to  an  autocorrelation  function  R(t),  is  given  by 


<^(f) 


R(x)e 


-j2irft 


dx. 


...  (20) 


It  is  assumed  that  for  the  discrete  system  being  considered,  the 
spectral  power  density  can  be  similarly  calculated  using  the  relationship 


(N-1) 

<t>(f)  = 2 SK^At 

n=-(N-l) 


...  (21) 


This  may  also  be  represented  as 


(N-1) 

♦ (f)  = 2 SK^At 

n=- (N-1) 


R (n,N)cos(2irfnAt)  , ...  (22) 


since  R (n,N)  is  an  even  function. 

For  the  actual  calculations  a modified  form  of  equation  (23)  was  used 

i.e . 


2M-2 


<^(f) 

2NSK^At 


1 

M-1 


2 Rl((m  - W-l),M)e-^2iTfmAt  ^ 


(23) 


m=l 


where  M = N + 1. 

The  lower  limit  of  m was  changed  from  0 to  1,  since  R^(n,M)  is  symmetrical 
about  R^(0,M),  the  number  of  coefficients  calculated  are  (2M-1).  For 
efficient  use  of  the  Fast  Fourier  Transform  (FFT)  the  R^(-(M-1),M)  term  is 
discarded  so  the  number  of  coefficients  is  a power  of  2.  This  was  found 
to  have  negligible  effect  on  the  results. 


Equation  (23)  has  been  evaluated  using  time  decimation  FFT  techniques 
for  various  values  of  N.  Sample  results  are  shovTi  in  the  log-log  plots  of 
Figures  2,  3,  4 and  5.  The  points  in  these  plots  appear  to  lie  within 
two  distinct  lines,  which  are  the  envelopes  of  the  oscillations  explained 


6 


in  Section  5 and  Illustrated  In  Figure  6.  The  continuous  straight  line 
shown  in  each  figure  is  the  least  squares  fitted  line  with  the  line  equation 
shown  on  the  figure.  The  frequency  exponent  a is  independent  of  At  but  is 
a function  of  N.  The  a values  calculated  using  the  least  squares  fits  are 
shown  in  Table  1. 


TABLE  1 

COMPARISON  OF  SPECTRAL  POWER  RESPONSE 
EXPONENTS  (f)~“  FOR  VARIATION  OF  N 


a 

256 

1.23 

512 

1.23 

1024 

1.22 

2048 

1.22 

The  fit  was  made  over  the  frequency  range  ^>5  f ^ t 'S 


It  is 


evident  that  a exceeds  unity,  but  decreases  as  N increases. 

Equation  (23)  can  be  used  directly  to  evaluate  the  low  frequency  (f-K)) 
power  spectrum  by  noting  that 


N-1 

E 

n=- (N-1) 


R (n,N)  « 4N  , 


...  (24) 


as  suggested  in  Appendix  II.  This  means  that 


lira 


»(f) 

2NSK^At 


= 4 


(25) 


This  low  frequency  asymptote  is  plotted  on  Figures  2,  3,  4,  5 and  6 as  well 
as  the  high  frequency  asymptote  for  the  continuous  case  derived  in 
Appendix  III. 
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As  f ->•  “ 


1 

2N(fAt) 


...  (26) 


»(f) 

2NSK^At 


5.  POWER  SPECTRUM  RIPPLE 


As  is  evident  in  Figure  6 and  suggested  in  Figures  2,  3,  4 and  5 by 
the  envelopes,  the  spectra  have  oscillatory  characteristics  superimposed  on 
a smoothly  varying  component.  It  is  suggested  by  Ternan  (6)  that  such 
behaviour  is  directly  attributable  to  the  truncation  of  the  impulse 
response.  In  Appendix  III  this  proposal  is  investigated  for  the  case  of  a 
continuous  system  with  an  abruptly  truncated  impulse  response.  It  is 
shovm  that  the  oscillatory  component  has  a repetition  "frequency"  in  the 
frequency  domain  of  NAt  and  that  the  ratio  of  oscillatory  amplitude  to 
basic  component  decreases  with  frequency  according  to  an  f“ ^ law.  This 
trend  is  also  expected  with  the  discrete  model  being  discussed. 

Figure  6 shows  an  evaluation  of  the  normalised  '-ersion  of  equation  (22). 


l.e. 


»(f) 

2NSK^At 


1 

N 


(N-1) 

E 

n=-  (N-1) 


R (n,N)cos (2irfnAt) 


(27) 


This  shows  the  fluctuations  in  the  low  frequency  end  of  the  spectrum,  and 
is  a direct  comparison  with  Figure  4. 

The  value  of  a = 1.21  shown  in  Figure  6 is  not  as  reliable  as  that  in 
Figure  4 because  of  the  small  number  of  points  used  in  the  least  squares 
fit . 


6.  CONCLUSION 


Spectral  analysis  of  the  process 


y^  = /2K(At) 


j-1 

E 

i=j-N 


X . 
1 


operating  on  samples  x.  from  a white  noise  generator  of  uniform  power 
spectral  density  S has  been  undertaken.  This  process  was  developed  to 
calculate  digital  values  of  quasi-flicker  noise  for  use  in  a computer 
simulation. 
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and  H(t)  = 2K(-2Trt)' 


-X 

X cos 


Now  It  may  be  shown,  using  standard  contour  integration  methods,  that 


-A 

X cos 


(■■¥) 


dx  = sin(^rA)  r(l-A)  , 


X cos 


dx  = o 


provided  o < X < 1. 

Substituting  the  expressions  from  (A7)  and  (A8)  into  (A5)  and  (A6) 
respectively,  gives 


H(t)  = 2K  sin(iTX)  r(l-X)(2irt)  , t>o 
H(t)  = o , t<o 


For  the  particular  condition  X = h,  and  noting  that  r(Js)  = /JT,  it 
follows  that 


H(t)  = /f  K t" 


...  (AlO) 


APPENDIX  II 


EVALUATION  OF  ER^(n.N) 


As  shown  In  the  text,  the  low  frequency  spectral  power  density  is  con- 
trolled by  the  value  of 

N-1 

A = ^ R^(n,N)  , 

n=- (N-1) 


where 


R^(n.N) 


, , -(N-l)<n^N-l 

(j+|n| ) 


No  standard  formula  appears  to  be  available  for  the  above  summation, 
although,  Macdonald  (8)  and  Harper  (9)  have  both  suggested  relationships  of 
the  form 


A-V4N  + 


C^(4N^  + 


N"*^) 


C2N 


-1 


+ C, 


) 


I 

I 

f 

i‘ 

f 


where  C2  and  are  constants,  for  evaluating  the  double  summation. 

However,  on  the  assumption  that  the  low  frequency  power  density  for  the 
discrete  case  and  the  continuous  case  as  analysed  in  Appendix  III  should  be 
closely  the  same,  then  it  is  expected  that  A''^4N  as  N-x».  This  proposition 
was  tested  numerically  for  various  N values  giving  the  results  below  : 


1 

N 

A 

4N 

10 

25.2 

40 

50 

162.6 

200 

100 

345.6 

400 

500 

1873.4 

2000 

1000 

3819.9 

4000 

It  is  concluded  that  the  proposition,  A = 4N,  is  accurate  within  5%  for 
N > 1000. 
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APPENDIX  III 


POWER  SPECTRUM  FOR  A CONTINUOUS  SYSTEM 


Consider  a linear  system  with  Impulse  response 


H(t)  = /2Kt 


-h  , 0 < t < T 


Its  transfer  function  is 


G(f) 


(u)e'^i27rfu  ^ 


which  on  substitution  of  (A6)  gives 


G(f)  = 


/2TrfT 


-h  -ju 
u e du 


LOW  FREQUENCY  BEHAVIOUR 

The  function  in  (A8)  can  be  written  in  series  form. 


G(f) 


E 


1 

(n  + h) 


(Z-rTfl)”"^  , 


from  which  it  follows  that 

G(o)  = 2/^  K 

The  output  spectral  power  density  for  a discrete  process  at  zero 
is  then 


<Ko)  = |g(o)|^  S 

= 8NAtK^  S , 


where  T,  the  duration  of  the  impulse  response,  is  replaced  by  NAt 
the  input  power  density. 


. . . (A6) 
...  (A7) 

...  (A8) 

...  (A9) 

. . . (AlO) 
frequency 

. ..  (All) 
and  S is 
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HIGH  FREQUENCY  BEHAVIOUR 


G(f),  from  (A8) , can  be  written 


G(f) 


K 


00  ^ 00 

-Ju  . K / -Jj  -ju  , 
u e du I u e du, 

/irf  Jl 

2 fT 


which,  using  results  from  Copson  (7),  gives 


G(f)  = — — [X  - JY] 


wf*^ 


(A12) 


where.  If  (fT)  Is  large  and  positive. 


X -\- 


(4^fT)3 


1.3.5  ^ 1.3. 5. 7. 9 ^ 


(AirfT)- 


and 


Y -V-  1 - -A-3  + lt3,_5..7 

(4TTfT)^  (4TrfT)^ 


(A13) 


If  we  consider  only 


2nfT  > 5 , 


then  Y ss  1 and  X « 0. 

The  output  spectral  power  density  4i(f)  Is  given  by 


^(f)  = SG(f)  G*(f)  . 

Substituting  (A12)  In  (A14)  with  the  restricted  values  of  X and 
normalised  power  spectrxim 


^(f) 

^ 1 

+ 1 

2NSK^At 

2N(fAt) 

4ir^N^(fAt)^ 

---I-  sin  - 2irN(fAt)l  , 

(fAt)  ' -I 

...  (A14) 
Y gives  the 


...  (A15) 
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where  as  before  NAt  = T. 


The  first  term  on  the  right  hand  side  represents  the  required  flicker 

-1  -2 
term  l.e.  f . The  second  term,  which  varies  with  f , and  the  third 

-3/2 

term,  which  oscillates  with  an  amplitude  proportional  to  f , can  be 
regarded  as  error  terms.  The  oscillatory  term  repeats  at  Intervals  of 
along  the  (fAt)  abscissa. 


2!|h- 
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FIG.  1 - Block  diagram  showing  the  Input  of  white  noise  x(t) 
to  a filter  with  impulse  response  H(t). 


DENS  I 


LOG  (FREQ  AT)  xiO 


OWER  DENS I 


-40,  -3S,  -30.  -25.  -20.  -1 

LOG  (FREQ.  AT)  xIO 

FIG.  4 


SPECTT^AL  POWER  DENSITY 


LOG  (FREQ  AT)  xIO 


LOG  (FREQ  AT)  X10 


(MRL-R-680) 


DISTRIBUTION  LIST 


MATERIALS  RESEARCH  LABORATORIES 


Chief  Superintendent 

Superintendent,  Physical  Chemistry  Division 

Dr.  C.I.  Sach 

Mr.  G.K.  Manzie 

Library 

Branch  Superintendent,  S.A.  Branch 

Librarian,  N.S.W.  Branch  (Through  Of f icer-in-Charge) 

Librarian,  S.A.  Branch 

Off icer-in-Charge,  Joint  Tropical  Research  Unit 

DEPARTMENT  OF  DEFENCE 

Chief  Defence  Scientist 

Executive  Controller,  ADSS 

Superintendent,  Central  Studies  Establishment 

Superintendent,  Defence  Science  Administration,  DSTO 

Superintendent,  Military  Advisers  Branch 

Head,  Laboratory  Programs  Branch 

Army  Scientific  Adviser 

Ik 

Air  Force  Scientific  Adviser 

Naval  Scientific  Adviser 

Chief  Superintendent,  Aeronautical  Research  Laboratories 

Director,  Weapons  Research  Establishment 

Senior  Librarian,  Weapons  Research  Establishment 

Librarian,  R.A.N.  Research  Laboratory 

Document  Exchange  Centre,  DLIS  (14  copies) 

1 

Principal  Librarian,  Campbell  Park  Library  ADSATIS  Annex 

■ 

Directorate  of  Quality  Assurance  (Air  Office) 

Head,  Engineering  Development  Establishment 

V • 

1 _ 

DEPARTMENT  OF  INDUSTRY  AND  COMMERCE 

NASA  Canberra  Office 
Head,  B.D.R.S.S.  (Aust.) 

OTHER  FEDERAL  AND  STATE  DEPARTMENTS  AND  INSTRUMENTALITIES 


The  Chief  Librarian,  Central  Library,  C.S.I.R.O. 

Australian  Atomic  Energy  Commission  Research  Establishment 


(MRL-R-680) 


DISTRIBUTION  LIST 
(Continued) 


MISCELLANEOUS  - OVERSEAS 


Defence  Scientific  & Technical  Representative,  Department  of 
Defence,  England 

Dr.  C.L.M.  Cottrell,  Assistant  Director/Armour  and  Materials, 
Military  Vehicles  and  Engineering  Establishment,  England 
Reports  Centre,  Directorate  of  Materials  Aviation,  England 
Library  - Exchange  Desk,  National  Bureau  of  Standards,  U.S.A. 
U.S.  Army  Standardization  Group,  Office  of  the  Scientific 
Standardization  Representative,  Canberra,  A.C.T. 

Senior  Standardization  Representative,  U.S.  Army  Standardization 
Group,  Canberra,  A.C.T. 

Chief,  Research  and  Development,  Defence  Scientific  Information 
Service,  Canada  (2  copies) 

The  Director,  Defence  Scientific  Information  and  Documentation 
Centre,  India 

Colonel  B.C.  Joshi,  Military,  Naval  and  Air  Adviser,  High 
Commission  of  India,  Red  Hill,  A.C.T. 

Director,  Defence  Research  Centre,  Malaysia 
Accessions  Department,  British  Library,  England 
Official  Publications  Section,  British  Library,  England 
Librarian,  Periodicals  Recording  Section,  National  Reference 
Library  of  Science  and  Invention,  England 
INSPEC:  Acquisition  Section,  Institution  of  Electrical 

Engineers,  England 

Overseas  Reports  Section,  Defence  Research  Information  Centre, 
England . 


I 


